# This script generates the analysis tables and figures

# Last run 19 December 2024

dat_conjoint <- read_csv("data/dat_conjoint.csv")

# Format characters as factors
dat_conjoint$outcome <- factor(dat_conjoint$outcome)
dat_conjoint$instrument <- factor(dat_conjoint$instrument)
dat_conjoint$coverage <- factor(dat_conjoint$coverage,levels = c("Narrow","Broad"))
dat_conjoint$revenue_use <- factor(dat_conjoint$revenue_use,levels = c("Green infrastructure","Compensate vulnerable","Compensate population"),
                                   labels = c("Green infrastructure","Compensate vulnerable","Compensate population"))

# 1. Construct conjoint design -----------------------------
attribute_list <- list()
attribute_list[["instrument"]] <- c("Carbon tax","Emissions trading")
attribute_list[["coverage"]] <- c("Narrow","Broad")
attribute_list[["revenue_use"]] <- c("Green infrastructure","Compensate vulnerable","Compensate population")

# Randomization constraints in the conjoint design (none)
constraint_list <- list()

# Make design
conjointdesign <- makeDesign(type = "constraints",attribute.levels = attribute_list,constraints = constraint_list)

# Default weights (uniform)

# 2. Main results -----------------------------

## 2.1 Table B1: Effective, Feasible, and Propose AMCEs -------------------------------------

# Run models
m_rq1 <- amce(selected_profile ~ instrument + coverage + revenue_use, 
              data = dat_conjoint %>% filter(outcome == "effective"), cluster = TRUE, respondent.id = "response_id",design = conjointdesign)
m_rq1_sum <- summary(m_rq1)

m_rq2 <- amce(selected_profile ~ instrument + coverage + revenue_use, 
              data = dat_conjoint %>% filter(outcome == "feasible"), cluster = TRUE, respondent.id = "response_id",design = conjointdesign)
m_rq2_sum <- summary(m_rq2)

m_rq3 <- amce(selected_profile ~ instrument + coverage + revenue_use, 
              data = dat_conjoint %>% filter(outcome == "propose"), cluster = TRUE, respondent.id = "response_id",design = conjointdesign)
m_rq3_sum <- summary(m_rq3)

mods_b1 <- list(m_rq1_sum,m_rq2_sum,m_rq3_sum)

# Generate table
tab_b1_vecs <- tab_func(mods_b1,amce_acie_resp = "amce")

tab_b1 <- tibble(attribute = c(rep("\\textbf{Coverage}",2),rep("\\textbf{Instrument}",2),rep("\\textbf{Revenue use}",4)),
                 vars = paste("\\multicolumn{1}{l}{",c(rep(mods_b1[[1]]$amce$Level,each = 2)),"}",sep = ""),
                 "eff"=NA_character_,"eff_adj"=NA_character_,"fea"=NA_character_,"fea_adj"=NA_character_,"pro"=NA_character_,"pro_adj"=NA_character_)
tab_b1[seq(1,nrow(tab_b1),by=2),-c(1:2,seq(4,ncol(tab_b1),by=2))] <- matrix(tab_b1_vecs$coefs,ncol = (ncol(tab_b1)-2)/2)
tab_b1[seq(2,nrow(tab_b1),by=2),-c(1:2,seq(4,ncol(tab_b1),by=2))] <- matrix(tab_b1_vecs$se,ncol = (ncol(tab_b1)-2)/2)
tab_b1[seq(1,nrow(tab_b1),by=2),-c(1:2,seq(3,ncol(tab_b1),by=2))] <- matrix(tab_b1_vecs$coefs_adj,ncol = (ncol(tab_b1)-2)/2)
tab_b1[seq(2,nrow(tab_b1),by=2),-c(1:2,seq(3,ncol(tab_b1),by=2))] <- matrix(tab_b1_vecs$se_adj,ncol = (ncol(tab_b1)-2)/2)
tab_b1[-c(1,3,5),1] <- ""
tab_b1[seq(2,nrow(tab_b1),by=2),2] <- ""
tab_b1 <- tab_b1[c(3:4,1:2,5:nrow(tab_b1)),]

Nobs <- c("\\midrule \\multicolumn{1}{c}{N(observations)}","",rep(m_rq1_sum$samplesize_estimates,2),rep(m_rq2_sum$samplesize_estimates,2),rep(m_rq3_sum$samplesize_estimates,2))
Nresp <- c("\\multicolumn{1}{c}{N(respondents)}","",rep(m_rq1_sum$respondents,2),rep(m_rq2_sum$respondents,2),rep(m_rq3_sum$respondents,2))
Ntab <- rbind(Nobs,Nresp); colnames(Ntab) <- colnames(tab_b1)
tab_b1 <- rbind(tab_b1,Ntab)

addtorow <- list(); addtorow$pos <- list(nrow(tab_b1)); addtorow$command <- c("")
ptable(tab_b1, file = "output/tb1_main.tex")

## 2.2 Figure 2: Effectiveness and Feasibility: Average Marginal Component Effects -------------------------------------

# Reformat table to plot data
plotdat_main <- tab_to_plotdat_func(tab_b1,adjustment = "holm")

plotdat_main <- plotdat_main %>%
  mutate(color = case_when(attribute == "Instrument" ~ "#EB657E",
                           attribute == "Coverage" ~ "#8F2A68",
                           attribute == "Revenue use" ~ "#3C2C49"),
         vars = glue("<span style='color:{color}'>{vars}"))

# Add base categories
base_categories <- tibble(attribute = rep(unique(plotdat_main$attribute),each = length(unique(plotdat_main$attribute))),
                          vars = rep(c("Carbon tax","Narrow","Green Infrastructure"), each = length(unique(plotdat_main$attribute))),
                          type = rep(c("Effective","Feasible","Propose"),length(unique(plotdat_main$attribute))),
                          est = 0,
                          se = 0,
                          color = rep(c("#EB657E","#8F2A68","#3C2C49"), each = length(unique(plotdat_main$attribute)))) %>%
  mutate(vars = glue("<span style='color:{color}'>**{attribute}**<br>{vars}</span>"))


# Merge
plotdat_main <- bind_rows(plotdat_main, base_categories)
levels_main <- factor(unique(plotdat_main$vars), levels = unique(plotdat_main$vars)[c(4,3,7,2,6,1,5)])

plotdat_main <- plotdat_main %>%
  mutate(attribute = factor(attribute, levels = unique(attribute)),
         vars = factor(vars,levels = levels(levels_main)),
         type = factor(type, levels = unique(type)),
         est = est/100,
         ll_95 = est - (qnorm(0.975)*se)/100,
         ul_95 = est + (qnorm(0.975)*se)/100,
         ll_90 = est - (qnorm(0.95)*se)/100,
         ul_90 = est + (qnorm(0.95)*se)/100,
  )

p_f2 <- ggplot(plotdat_main %>% filter(type %in% c("Effective","Feasible"))) +
  geom_hline(yintercept = 0,color = "grey") +
  geom_point(aes(x = vars, y = est, color = attribute, group = attribute)) +
  geom_errorbar(aes(x = vars, y = est, ymin = ll_90, ymax = ul_90,
                    color = attribute, group = attribute), width = 0, linewidth = 1) +
  geom_errorbar(aes(x = vars, y = est, ymin = ll_95, ymax = ul_95,
                    color = attribute, group = attribute), width = 0) +
  scale_y_continuous(labels = percent_format()) +
  scale_color_manual(values = unique(plotdat_main$color)) +
  coord_flip() +
  facet_wrap(~type,ncol = 3) +
  labs(x = NULL, y = "Change in probability of selection") +
  theme_bw() +
  theme(
    axis.text.y = element_markdown(size = 10),
    axis.text.x = element_text(size = 6),
    axis.title.x = element_text(face = "bold",size = 10),
    plot.caption = element_markdown(lineheight = 1.2),panel.grid = element_blank(),strip.text = element_text(face = "bold"),
    legend.position = "none",panel.spacing = unit(0.5,units = "cm"),
  )

ggsave(p_f2,file = "output/p_f2.png",width = 6, height = 4,units = "in")

## 2.3 Table D1: Association between effectiveness, feasibility, and the probability of policy proposal -------------------------------------
# Equivalent to lm_robust(propose ~ effective * feasible,clusters = response_contest)

tab_td1 <- dat_conjoint %>%
  pivot_wider(id_cols = c("response_id","contest_no","profile"),names_from = outcome,values_from = selected_profile) %>%
  mutate(response_contest = paste0(response_id, "_",contest_no)) %>%
  na.omit() %>%
  group_by(effective, feasible) %>%
  summarize(p = paste0(round(sum(propose)/n()*100,digits = 1),"\\%"),.groups = 'drop') %>%
  pivot_wider(names_from = "feasible",values_from = "p") %>%
  rename(`Probability of Proposal` = effective, `Less Feasible` = `0`, `More Feasible` = `1`) %>%
  mutate(`Probability of Proposal` = c("Less Effective","More Effective"))

addtorow <- list(); addtorow$pos <- list(2); addtorow$command <- c("")
ptable(tab_td1, file = "output/td1_crosstab.tex")

## 2.4 Figure 3: Proposed Choice: Average Marginal Component Effects -------------------------------------

p_f3 <- ggplot(plotdat_main %>% filter(type %in% c("Propose"))) +
  geom_hline(yintercept = 0,color = "grey") +
  geom_point(aes(x = vars, y = est, color = attribute, group = attribute)) +
  geom_errorbar(aes(x = vars, y = est, ymin = ll_90, ymax = ul_90,
                    color = attribute, group = attribute), width = 0, linewidth = 1) +
  geom_errorbar(aes(x = vars, y = est, ymin = ll_95, ymax = ul_95,
                    color = attribute, group = attribute), width = 0) +
  scale_y_continuous(labels = percent_format(),breaks = seq(-.3,.2,by=.1),limits = c(-.3,.2)) +
  scale_color_manual(values = unique(plotdat_main$color)) +
  coord_flip() +
  facet_wrap(~type,ncol = 3) +
  labs(x = NULL, y = "Change in probability of selection") +
  theme_bw() +
  theme(
    axis.text.y = element_markdown(size = 10),
    axis.text.x = element_text(size = 6),
    axis.title.x = element_text(face = "bold",size = 10),
    plot.caption = element_markdown(lineheight = 1.2),panel.grid = element_blank(),strip.text = element_text(face = "bold"),
    legend.position = "none",panel.spacing = unit(0.5,units = "cm"),
  )

ggsave(p_f3,file = "output/p_f3.png",width = 4, height = 4,units = "in")

# 3. Appendix C: Interactions between policy design attributes -------------------------------------

## Instrument x Revenue use
m_rq5_eff_ixr <- amce(selected_profile ~ instrument + coverage + revenue_use + instrument:revenue_use, 
                      data = dat_conjoint %>% filter(outcome == "effective"), cluster = TRUE, respondent.id = "response_id",design = conjointdesign)
m_rq5_eff_ixr_sum <- summary(m_rq5_eff_ixr)
m_rq5_fea_ixr <- amce(selected_profile ~ instrument + coverage + revenue_use + instrument:revenue_use, 
                      data = dat_conjoint %>% filter(outcome == "feasible"), cluster = TRUE, respondent.id = "response_id",design = conjointdesign)
m_rq5_fea_ixr_sum <- summary(m_rq5_fea_ixr)
m_rq5_pro_ixr <- amce(selected_profile ~ instrument + coverage + revenue_use + instrument:revenue_use, 
                      data = dat_conjoint %>% filter(outcome == "propose"), cluster = TRUE, respondent.id = "response_id",design = conjointdesign)
m_rq5_pro_ixr_sum <- summary(m_rq5_pro_ixr)

## Coverage x Revenue use
m_rq5_eff_cxr <- amce(selected_profile ~ instrument + coverage + revenue_use + coverage:revenue_use, 
                      data = dat_conjoint %>% filter(outcome == "effective"), cluster = TRUE, respondent.id = "response_id",design = conjointdesign)
m_rq5_eff_cxr_sum <- summary(m_rq5_eff_cxr)
m_rq5_fea_cxr <- amce(selected_profile ~ instrument + coverage + revenue_use + coverage:revenue_use, 
                      data = dat_conjoint %>% filter(outcome == "feasible"), cluster = TRUE, respondent.id = "response_id",design = conjointdesign)
m_rq5_fea_cxr_sum <- summary(m_rq5_fea_cxr)
m_rq5_pro_cxr <- amce(selected_profile ~ instrument + coverage + revenue_use + coverage:revenue_use, 
                      data = dat_conjoint %>% filter(outcome == "propose"), cluster = TRUE, respondent.id = "response_id",design = conjointdesign)
m_rq5_pro_cxr_sum <- summary(m_rq5_pro_cxr)

## Instrument x Coverage
m_rq5_eff_ixc <- amce(selected_profile ~ instrument + coverage + revenue_use + instrument:coverage, 
                      data = dat_conjoint %>% filter(outcome == "effective"), cluster = TRUE, respondent.id = "response_id",design = conjointdesign)
m_rq5_eff_ixc_sum <- summary(m_rq5_eff_ixc)
m_rq5_fea_ixc <- amce(selected_profile ~ instrument + coverage + revenue_use + instrument:coverage, 
                      data = dat_conjoint %>% filter(outcome == "feasible"), cluster = TRUE, respondent.id = "response_id",design = conjointdesign)
m_rq5_fea_ixc_sum <- summary(m_rq5_fea_ixc)
m_rq5_pro_ixc <- amce(selected_profile ~ instrument + coverage + revenue_use + instrument:coverage, 
                      data = dat_conjoint %>% filter(outcome == "propose"), cluster = TRUE, respondent.id = "response_id",design = conjointdesign)
m_rq5_pro_ixc_sum <- summary(m_rq5_pro_ixc)

## 3.1 Table C1: ACIEs by policy design attribute -------------------------------------
mods_tc1 <- list(m_rq5_eff_cxr_sum,m_rq5_eff_ixc_sum,m_rq5_eff_ixr_sum,
                 m_rq5_fea_cxr_sum,m_rq5_fea_ixc_sum,m_rq5_fea_ixr_sum,
                 m_rq5_pro_cxr_sum,m_rq5_pro_ixc_sum,m_rq5_pro_ixr_sum)

tab_tc1_vecs <- tab_func(mods_tc1,amce_acie_resp = "acie")
tab_tc1_vecs$se_adj <- if_else(tab_tc1_vecs$se_adj=="(Inf)","\\multicolumn{1}{c}{\\textit{N.S.}}",tab_tc1_vecs$se_adj)

tab_tc1 <- tibble(vars = paste("\\multicolumn{1}{l}{",c(rep(as.vector(unlist(sapply(mods_tc1[1:3],function(x){x$acie$Level}))),each = 2)),"}",sep = ""),
                  "eff"=NA_character_,"eff_adj"=NA_character_,"fea"=NA_character_,"fea_adj"=NA_character_,"pro"=NA_character_,"pro_adj"=NA_character_)
tab_tc1[seq(1,nrow(tab_tc1),by=2),-c(1,seq(3,ncol(tab_tc1),by=2))] <- matrix(tab_tc1_vecs$coefs,ncol = (ncol(tab_tc1)-1)/2)
tab_tc1[seq(2,nrow(tab_tc1),by=2),-c(1,seq(3,ncol(tab_tc1),by=2))] <- matrix(tab_tc1_vecs$se,ncol = (ncol(tab_tc1)-1)/2)
tab_tc1[seq(1,nrow(tab_tc1),by=2),-c(1,seq(2,ncol(tab_tc1),by=2))] <- matrix(tab_tc1_vecs$coefs_adj,ncol = (ncol(tab_tc1)-1)/2)
tab_tc1[seq(2,nrow(tab_tc1),by=2),-c(1,seq(2,ncol(tab_tc1),by=2))] <- matrix(tab_tc1_vecs$se_adj,ncol = (ncol(tab_tc1)-1)/2)
tab_tc1[seq(2,nrow(tab_tc1),by=2),1] <- ""
tab_tc1$vars <- str_replace(tab_tc1$vars,pattern = "Broad[:]Emissions trading",replacement = "Emissions trading:Broad")
tab_tc1 <- tab_tc1[c(5:10,1:4),]

Nobs <- c("\\midrule \\multicolumn{1}{l}{N(observations)}",rep(m_rq5_eff_cxr_sum$samplesize_estimates,2),rep(m_rq5_fea_cxr_sum$samplesize_estimates,2),rep(m_rq5_pro_cxr_sum$samplesize_estimates,2))
Nresp <- c("\\multicolumn{1}{l}{N(respondents)}",rep(m_rq5_eff_cxr_sum$respondents,2),rep(m_rq5_fea_cxr_sum$respondents,2),rep(m_rq5_pro_cxr_sum$respondents,2))
Ntab <- rbind(Nobs,Nresp); colnames(Ntab) <- colnames(tab_tc1)
tab_tc1 <- rbind(tab_tc1,Ntab)

addtorow <- list(); addtorow$pos <- list(nrow(tab_tc1)); addtorow$command <- c("")
ptable(tab_tc1, file = "output/tc1_interactions.tex")

# 4. Appendix E: Heterogeneous treatment effects -------------------------------------

## 4.1 Iteratively drop each country  -------------------------------------
### 4.1.1 Figure E1: Country leave-one-out analysis -------------------------------------
result_loo <- tibble(type = NA_character_, country_out = NA_character_, attribute = NA_character_, vars = NA_character_, est = NA_real_, se = NA_real_, se_corrected = NA_real_,.rows = 0)

for(c in sort(unique(dat_conjoint$country))){
  dat_conjoint_loo <- dat_conjoint %>%
    filter(country %in% c == F)
  
  m_eff_loo <- amce(selected_profile ~ instrument + coverage + revenue_use, 
                    data = dat_conjoint_loo %>% filter(outcome == "effective"), cluster = TRUE, respondent.id = "response_id",design = conjointdesign)
  result_loo <- bind_rows(result_loo,
                          tibble(type = "effective", country_out = c, 
                                 attribute = summary(m_eff_loo)$amce$Attribute, vars = summary(m_eff_loo)$amce$Level,
                                 est = summary(m_eff_loo)$amce$Estimate, se = summary(m_eff_loo)$amce$`Std. Err`,
                                 se_corrected = abs(summary(m_eff_loo)$amce$Estimate/qnorm(1-(p.adjust(p = summary(m_eff_loo)$amce$`Pr(>|z|)`,method = "holm")/2)))))
  
  m_fea_loo <-  amce(selected_profile ~ instrument + coverage + revenue_use, 
                     data = dat_conjoint_loo %>% filter(outcome == "feasible"), cluster = TRUE, respondent.id = "response_id",design = conjointdesign)
  result_loo <- bind_rows(result_loo,
                          tibble(type = "feasible", country_out = c, 
                                 attribute = summary(m_fea_loo)$amce$Attribute, vars = summary(m_fea_loo)$amce$Level,
                                 est = summary(m_fea_loo)$amce$Estimate, se = summary(m_fea_loo)$amce$`Std. Err`,
                                 se_corrected = abs(summary(m_fea_loo)$amce$Estimate/qnorm(1-(p.adjust(p = summary(m_fea_loo)$amce$`Pr(>|z|)`,method = "holm")/2)))))
  
  m_pro_loo <- amce(selected_profile ~ instrument + coverage + revenue_use, 
                    data = dat_conjoint_loo %>% filter(outcome == "propose"), cluster = TRUE, respondent.id = "response_id",design = conjointdesign)
  result_loo <- bind_rows(result_loo,
                          tibble(type = "propose", country_out = c, 
                                 attribute = summary(m_pro_loo)$amce$Attribute, vars = summary(m_pro_loo)$amce$Level,
                                 est = summary(m_pro_loo)$amce$Estimate, se = summary(m_pro_loo)$amce$`Std. Err`,
                                 se_corrected = abs(summary(m_pro_loo)$amce$Estimate/qnorm(1-(p.adjust(p = summary(m_pro_loo)$amce$`Pr(>|z|)`,method = "holm")/2)))))
  
}

# Figure E1: Country leave-one-out AMCEs
plotdat_loo <- result_loo %>%
  mutate(type = str_to_title(type),
         attribute = str_replace(str_to_title(attribute),pattern = "_",replacement = " "),
         color = case_when(attribute == "Instrument" ~ "#EB657E",
                           attribute == "Coverage" ~ "#8F2A68",
                           attribute == "Revenue use" ~ "#3C2C49"),
         vars = glue("<span style='color:{color}'>{vars}"))

levels_loo <- factor(unique(plotdat_loo$vars), levels = unique(plotdat_loo$vars)[c(4,3,7,2,6,1,5)])

plotdat_loo <- plotdat_loo %>%
  mutate(attribute = factor(attribute, levels = unique(attribute)),
         vars = factor(vars,levels = levels(levels_loo)),
         type = factor(type, levels = unique(type)),
         est = est,
         ll_95 = est - (qnorm(0.975)*se),
         ul_95 = est + (qnorm(0.975)*se),
         ll_90 = est - (qnorm(0.95)*se),
         ul_90 = est + (qnorm(0.95)*se),
  )

p_fe1_loo <- ggplot(plotdat_loo) +
  geom_hline(yintercept = 0,color = "grey") +
  geom_point(data = plotdat_main, aes(x = vars, y = est, color = attribute, group = attribute)) +
  geom_point(aes(x = vars, y = est, color = attribute, group = attribute),position = position_jitter(width = 0.1,seed = 12123), alpha = 0.05) +
  geom_errorbar(data = plotdat_main, aes(x = vars, y = est, ymin = ll_90, ymax = ul_90, 
                                       color = attribute, group = attribute), width = 0, linewidth = 1) +
  geom_errorbar(aes(x = vars, y = est, ymin = ll_90, ymax = ul_90, 
                    color = attribute, group = attribute), width = 0, linewidth = 1,position = position_jitter(width = 0.1,seed = 12123), alpha = 0.05) +
  geom_errorbar(data = plotdat_main, aes(x = vars, y = est, ymin = ll_95, ymax = ul_95, 
                                       color = attribute, group = attribute), width = 0) +
  geom_errorbar(aes(x = vars, y = est, ymin = ll_95, ymax = ul_95, 
                    color = attribute, group = attribute), width = 0,position = position_jitter(width = 0.1,seed = 12123), alpha = 0.05) +
  scale_y_continuous(labels = percent_format(),breaks = seq(-.3,.2,by=.1),limits = c(-.3,.2)) +
  scale_color_manual(values = unique(plotdat_main$color)) +
  coord_flip() +
  facet_wrap(~type,ncol = 3) +
  labs(x = NULL, y = "Change in probability of selection") +
  theme_bw() +
  theme(
    axis.text.y = element_markdown(size = 10),
    axis.text.x = element_text(size = 6),
    axis.title.x = element_text(face = "bold",size = 10),
    plot.caption = element_markdown(lineheight = 1.2),panel.grid = element_blank(),strip.text = element_text(face = "bold"),
    legend.position = "none",panel.spacing = unit(0.25,units = "cm"),
  )

ggsave(p_fe1_loo, file = "output/p_e1_loo.png",width = 8, height = 4,units = "in")

### 4.1.2 Figure E2: Country leave-two-out analysis -------------------------------------
result_lto <- tibble(type = NA_character_, country_out_1 = NA_character_, country_out_2 = NA_character_, attribute = NA_character_, vars = NA_character_, est = NA_real_, se = NA_real_, se_corrected = NA_real_,.rows = 0)

for(c1 in sort(unique(dat_conjoint$country))){
  for(c2 in sort(unique(dat_conjoint$country[-which(dat_conjoint$country==c1)]))){
    dat_conjoint_lto <- dat_conjoint %>%
      filter(country != c1 & country != c2)
    
    m_eff_lto <- amce(selected_profile ~ instrument + coverage + revenue_use, 
                      data = dat_conjoint_lto %>% filter(outcome == "effective"), cluster = TRUE, respondent.id = "response_id",design = conjointdesign)
    result_lto <- bind_rows(result_lto,
                            tibble(type = "effective", country_out_1 = c1, country_out_2 = c2, 
                                   attribute = summary(m_eff_lto)$amce$Attribute, vars = summary(m_eff_lto)$amce$Level,
                                   est = summary(m_eff_lto)$amce$Estimate, se = summary(m_eff_lto)$amce$`Std. Err`,
                                   se_corrected = abs(summary(m_eff_lto)$amce$Estimate/qnorm(1-(p.adjust(p = summary(m_eff_lto)$amce$`Pr(>|z|)`,method = "holm")/2)))))
    
    m_fea_lto <-  amce(selected_profile ~ instrument + coverage + revenue_use, 
                       data = dat_conjoint_lto %>% filter(outcome == "feasible"), cluster = TRUE, respondent.id = "response_id",design = conjointdesign)
    result_lto <- bind_rows(result_lto,
                            tibble(type = "feasible", country_out_1 = c1, country_out_2 = c2, 
                                   attribute = summary(m_fea_lto)$amce$Attribute, vars = summary(m_fea_lto)$amce$Level,
                                   est = summary(m_fea_lto)$amce$Estimate, se = summary(m_fea_lto)$amce$`Std. Err`,
                                   se_corrected = abs(summary(m_fea_lto)$amce$Estimate/qnorm(1-(p.adjust(p = summary(m_fea_lto)$amce$`Pr(>|z|)`,method = "holm")/2)))))
    
    m_pro_lto <- amce(selected_profile ~ instrument + coverage + revenue_use, 
                      data = dat_conjoint_lto %>% filter(outcome == "propose"), cluster = TRUE, respondent.id = "response_id",design = conjointdesign)
    result_lto <- bind_rows(result_lto,
                            tibble(type = "propose", country_out_1 = c1, country_out_2 = c2, 
                                   attribute = summary(m_pro_lto)$amce$Attribute, vars = summary(m_pro_lto)$amce$Level,
                                   est = summary(m_pro_lto)$amce$Estimate, se = summary(m_pro_lto)$amce$`Std. Err`,
                                   se_corrected = abs(summary(m_pro_lto)$amce$Estimate/qnorm(1-(p.adjust(p = summary(m_pro_lto)$amce$`Pr(>|z|)`,method = "holm")/2)))))
  }
}

# Remove rows in which the pairing (country_out_1, country_out_2) is duplicated as (country_out_2, country_out_1)
result_lto <- result_lto %>% filter(country_out_1 < country_out_2)

# Figure E2: Country leave-two-out AMCEs
plotdat_lto <- result_lto %>%
  mutate(type = str_to_title(type),
         attribute = str_replace(str_to_title(attribute),pattern = "_",replacement = " "),
         color = case_when(attribute == "Instrument" ~ "#EB657E",
                           attribute == "Coverage" ~ "#8F2A68",
                           attribute == "Revenue use" ~ "#3C2C49"),
         vars = glue("<span style='color:{color}'>{vars}"))

levels_lto <- factor(unique(plotdat_lto$vars), levels = unique(plotdat_lto$vars)[c(4,3,7,2,6,1,5)])

plotdat_lto <- plotdat_lto %>%
  mutate(attribute = factor(attribute, levels = unique(attribute)),
         vars = factor(vars,levels = levels(levels_lto)),
         type = factor(type, levels = unique(type)),
         est = est,
         ll_95 = est - (qnorm(0.975)*se),
         ul_95 = est + (qnorm(0.975)*se),
         ll_90 = est - (qnorm(0.95)*se),
         ul_90 = est + (qnorm(0.95)*se),
  )

p_fe2_lto <- ggplot(plotdat_lto) +
  geom_hline(yintercept = 0,color = "grey") +
  geom_point(data = plotdat_main, aes(x = vars, y = est, color = attribute, group = attribute)) +
  geom_point(aes(x = vars, y = est, color = attribute, group = attribute),position = position_jitter(width = 0.1,seed = 12123), alpha = 0.02) +
  geom_errorbar(data = plotdat_main, aes(x = vars, y = est, ymin = ll_90, ymax = ul_90, 
                                         color = attribute, group = attribute), width = 0, linewidth = 1) +
  geom_errorbar(aes(x = vars, y = est, ymin = ll_90, ymax = ul_90, 
                    color = attribute, group = attribute), width = 0, linewidth = 1,position = position_jitter(width = 0.1,seed = 12123), alpha = 0.02) +
  geom_errorbar(data = plotdat_main, aes(x = vars, y = est, ymin = ll_95, ymax = ul_95, 
                                         color = attribute, group = attribute), width = 0) +
  geom_errorbar(aes(x = vars, y = est, ymin = ll_95, ymax = ul_95, 
                    color = attribute, group = attribute), width = 0,position = position_jitter(width = 0.1,seed = 12123), alpha = 0.02) +
  scale_y_continuous(labels = percent_format(),breaks = seq(-.35,.2,by=.1),limits = c(-.35,.2)) +
  scale_color_manual(values = unique(plotdat_main$color)) +
  coord_flip() +
  facet_wrap(~type,ncol = 3) +
  labs(x = NULL, y = "Change in probability of selection") +
  theme_bw() +
  theme(
    axis.text.y = element_markdown(size = 10),
    axis.text.x = element_text(size = 6),
    axis.title.x = element_text(face = "bold",size = 10),
    plot.caption = element_markdown(lineheight = 1.2),panel.grid = element_blank(),strip.text = element_text(face = "bold"),
    legend.position = "none",panel.spacing = unit(0.25,units = "cm"),
  )

ggsave(p_fe2_lto, file = "output/p_e2_lto.png",width = 8, height = 4,units = "in")

## 4.2 Table E1: Re-weighting by sample frame --------------
# Load sample frame target proportions
sample_frame_fractions <- read_csv("data/sample_frame_fractions.csv")

targets <- with(sample_frame_fractions,list(
  gender = wpct(gender, fraction),
  in_government = wpct(ingovernment, fraction)
))

dat_conjoint_weights <- data.frame(dat_conjoint %>%
  select(response_id, gender, in_government) %>%
  na.omit() %>%
  distinct() %>%
  mutate(gender = factor(gender),
         in_government = factor(in_government)))

# Use automated raking, iterating through cases and adjusting weight until sample and population distributions are aligned along variables of interest
myweights <- anesrake(targets,dat_conjoint_weights, caseid = dat_conjoint_weights$response_id,type = "nolim")

myweights <- left_join(dat_conjoint %>% select(response_id), 
                       tibble(response_id = myweights$caseid,weight = myweights$weightvec),
                       by = join_by(response_id)) %>%
  distinct()

dat_conjoint <- left_join(dat_conjoint,myweights,by = join_by(response_id))

m_rq1_te1 <- amce(selected_profile ~ instrument + coverage + revenue_use, 
              data = dat_conjoint %>% select(response_id,outcome,selected_profile, instrument,coverage,revenue_use,weight) %>% na.omit() %>% filter(outcome == "effective"), 
              cluster = TRUE, respondent.id = "response_id",design = conjointdesign,weights = "weight")
m_rq1_te1_sum <- summary(m_rq1_te1)

m_rq2_te1 <- amce(selected_profile ~ instrument + coverage + revenue_use, 
                data = dat_conjoint %>% select(response_id,outcome,selected_profile, instrument,coverage,revenue_use,weight) %>% na.omit() %>% filter(outcome == "feasible"), 
                cluster = TRUE, respondent.id = "response_id",design = conjointdesign,weights = "weight")
m_rq2_te1_sum <- summary(m_rq2_te1)

m_rq3_te1 <- amce(selected_profile ~ instrument + coverage + revenue_use, 
                data = dat_conjoint %>% select(response_id,outcome,selected_profile, instrument,coverage,revenue_use,weight) %>% na.omit() %>% filter(outcome == "propose"),
                cluster = TRUE, respondent.id = "response_id",design = conjointdesign,weights = "weight")
m_rq3_te1_sum <- summary(m_rq3_te1)

mods_t_te1 <- list(m_rq1_te1_sum,m_rq2_te1_sum,m_rq3_te1_sum)

tab_t_te1_vecs <- tab_func(mods_t_te1,amce_acie_resp = "amce")

tab_t_te1 <- tibble(attribute = c(rep("\\textbf{Coverage}",2),rep("\\textbf{Instrument}",2),rep("\\textbf{Revenue use}",4)),
                  vars = paste("\\multicolumn{1}{l}{",c(rep(mods_t_te1[[1]]$amce$Level,each = 2)),"}",sep = ""),
                  "eff"=NA_character_,"eff_adj"=NA_character_,"fea"=NA_character_,"fea_adj"=NA_character_,"pro"=NA_character_,"pro_adj"=NA_character_)
tab_t_te1[seq(1,nrow(tab_t_te1),by=2),-c(1:2,seq(4,ncol(tab_t_te1),by=2))] <- matrix(tab_t_te1_vecs$coefs,ncol = (ncol(tab_t_te1)-2)/2)
tab_t_te1[seq(2,nrow(tab_t_te1),by=2),-c(1:2,seq(4,ncol(tab_t_te1),by=2))] <- matrix(tab_t_te1_vecs$se,ncol = (ncol(tab_t_te1)-2)/2)
tab_t_te1[seq(1,nrow(tab_t_te1),by=2),-c(1:2,seq(3,ncol(tab_t_te1),by=2))] <- matrix(tab_t_te1_vecs$coefs_adj,ncol = (ncol(tab_t_te1)-2)/2)
tab_t_te1[seq(2,nrow(tab_t_te1),by=2),-c(1:2,seq(3,ncol(tab_t_te1),by=2))] <- matrix(tab_t_te1_vecs$se_adj,ncol = (ncol(tab_t_te1)-2)/2)
tab_t_te1[-c(1,3,5),1] <- ""
tab_t_te1[seq(2,nrow(tab_t_te1),by=2),2] <- ""
tab_t_te1 <- tab_t_te1[c(3:4,1:2,5:nrow(tab_t_te1)),]

Nobs <- c("\\midrule \\multicolumn{1}{c}{N(observations)}","",rep(m_rq1_te1_sum$samplesize_estimates,2),rep(m_rq2_te1_sum$samplesize_estimates,2),rep(m_rq3_te1_sum$samplesize_estimates,2))
Nresp <- c("\\multicolumn{1}{c}{N(respondents)}","",rep(m_rq1_te1_sum$respondents,2),rep(m_rq2_te1_sum$respondents,2),rep(m_rq3_te1_sum$respondents,2))
Ntab <- rbind(Nobs,Nresp); colnames(Ntab) <- colnames(tab_t_te1)
tab_t_te1 <- rbind(tab_t_te1,Ntab)

addtorow <- list(); addtorow$pos <- list(nrow(tab_t_te1)); addtorow$command <- c("")
ptable(tab_t_te1, "output/te1_weight.tex")


## 4.2 Table E2: Restricted sample to complete-responders -------------------------------------
all_qs <- dat_conjoint %>%
  group_by(response_id) %>%
  summarize(all_responded = if_else(sum(selected_profile)==18,1,0)) %>%
  filter(all_responded==1)

m_rq1_te2 <- amce(selected_profile ~ instrument + coverage + revenue_use, 
              data = dat_conjoint %>% filter(outcome == "effective") %>%
                filter(response_id %in% all_qs$response_id), cluster = TRUE, respondent.id = "response_id",design = conjointdesign)
m_rq1_sum_te2 <- summary(m_rq1_te2)

m_rq2_te2 <- amce(selected_profile ~ instrument + coverage + revenue_use, 
              data = dat_conjoint %>% filter(outcome == "feasible") %>%
                filter(response_id %in% all_qs$response_id), cluster = TRUE, respondent.id = "response_id",design = conjointdesign)
m_rq2_sum_te2 <- summary(m_rq2_te2)

m_rq3_te2 <- amce(selected_profile ~ instrument + coverage + revenue_use, 
              data = dat_conjoint %>% filter(outcome == "propose") %>%
                filter(response_id %in% all_qs$response_id), cluster = TRUE, respondent.id = "response_id",design = conjointdesign)
m_rq3_sum_te2 <- summary(m_rq3_te2)

mods_te2 <- list(m_rq1_sum_te2,m_rq2_sum_te2,m_rq3_sum_te2)

# Generate table
tab_te2_vecs <- tab_func(mods_te2,amce_acie_resp = "amce")

tab_te2 <- tibble(attribute = c(rep("\\textbf{Coverage}",2),rep("\\textbf{Instrument}",2),rep("\\textbf{Revenue use}",4)),
                  vars = paste("\\multicolumn{1}{l}{",c(rep(mods_te2[[1]]$amce$Level,each = 2)),"}",sep = ""),
                  "eff"=NA_character_,"eff_adj"=NA_character_,"fea"=NA_character_,"fea_adj"=NA_character_,"pro"=NA_character_,"pro_adj"=NA_character_)
tab_te2[seq(1,nrow(tab_te2),by=2),-c(1:2,seq(4,ncol(tab_te2),by=2))] <- matrix(tab_te2_vecs$coefs,ncol = (ncol(tab_te2)-2)/2)
tab_te2[seq(2,nrow(tab_te2),by=2),-c(1:2,seq(4,ncol(tab_te2),by=2))] <- matrix(tab_te2_vecs$se,ncol = (ncol(tab_te2)-2)/2)
tab_te2[seq(1,nrow(tab_te2),by=2),-c(1:2,seq(3,ncol(tab_te2),by=2))] <- matrix(tab_te2_vecs$coefs_adj,ncol = (ncol(tab_te2)-2)/2)
tab_te2[seq(2,nrow(tab_te2),by=2),-c(1:2,seq(3,ncol(tab_te2),by=2))] <- matrix(tab_te2_vecs$se_adj,ncol = (ncol(tab_te2)-2)/2)
tab_te2[-c(1,3,5),1] <- ""
tab_te2[seq(2,nrow(tab_te2),by=2),2] <- ""
tab_te2 <- tab_te2[c(3:4,1:2,5:nrow(tab_te2)),]

Nobs <- c("\\midrule \\multicolumn{1}{c}{N(observations)}","",rep(m_rq1_sum_te2$samplesize_estimates,2),rep(m_rq2_sum_te2$samplesize_estimates,2),rep(m_rq3_sum_te2$samplesize_estimates,2))
Nresp <- c("\\multicolumn{1}{c}{N(respondents)}","",rep(m_rq1_sum_te2$respondents,2),rep(m_rq2_sum_te2$respondents,2),rep(m_rq3_sum_te2$respondents,2))
Ntab <- rbind(Nobs,Nresp); colnames(Ntab) <- colnames(tab_te2)
tab_te2 <- rbind(tab_te2,Ntab)

addtorow <- list(); addtorow$pos <- list(nrow(tab_te2)); addtorow$command <- c("")
ptable(tab_te2, file = "output/te2_allqs.tex")


## 4.3 Carbon pricing context -------------------------------------

### 4.3.1 Table E3: Carbon pricing status -------------------------------------
m_hte8_eff <- amce(selected_profile ~ instrument + coverage + revenue_use + cp_status_factor + 
                     cp_status_factor:instrument + cp_status_factor:coverage + cp_status_factor:revenue_use, 
                   data = dat_conjoint %>% filter(outcome == "effective") %>% 
                     mutate(cp_status_factor = factor(cp_status_factor,levels = c(0,1,2),labels = c("None","Under consideration","Adopted"))) %>% 
                     select(response_id, selected_profile, instrument, coverage, revenue_use, cp_status_factor) %>% na.omit(),
                   cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                   respondent.varying = "cp_status_factor")
m_hte8_eff_sum <- summary(m_hte8_eff)

m_hte8_fea <- amce(selected_profile ~ instrument + coverage + revenue_use + cp_status_factor + 
                     cp_status_factor:instrument + cp_status_factor:coverage + cp_status_factor:revenue_use, 
                   data = dat_conjoint %>% filter(outcome == "feasible") %>% mutate(cp_status_factor = factor(cp_status_factor,levels = c(0,1,2),labels = c("None","Under consideration","Adopted"))) %>% 
                     select(response_id, selected_profile, instrument, coverage, revenue_use, cp_status_factor) %>% na.omit(),
                   cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                   respondent.varying = "cp_status_factor")
m_hte8_fea_sum <- summary(m_hte8_fea)

m_hte8_pro <- amce(selected_profile ~ instrument + coverage + revenue_use + cp_status_factor + 
                     cp_status_factor:instrument + cp_status_factor:coverage + cp_status_factor:revenue_use, 
                   data = dat_conjoint %>% filter(outcome == "propose") %>%
                     mutate(cp_status_factor = factor(cp_status_factor,levels = c(0,1,2),labels = c("None","Under consideration","Adopted"))) %>% 
                     select(response_id, selected_profile, instrument, coverage, revenue_use, cp_status_factor) %>% na.omit(),
                   cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                   respondent.varying = "cp_status_factor")
m_hte8_pro_sum <- summary(m_hte8_pro)

# Make table
mods_te3_hte8 <- list(m_hte8_eff_sum,m_hte8_fea_sum,m_hte8_pro_sum)
mods_te3_hte8_vecs <- tab_func(mods_te3_hte8,amce_acie_resp = "resp")
mods_te3_hte8_vecs$se_adj <- if_else(mods_te3_hte8_vecs$se_adj=="(Inf)","\\multicolumn{1}{c}{\\textit{N.S.}}",mods_te3_hte8_vecs$se_adj)
level_values_hte8 <- c("None","Under consideration","Adopted")

tab_te3_hte8 <- tab_hte_func(mods = mods_te3_hte8, mods_vecs = mods_te3_hte8_vecs, level_values = level_values_hte8)

# Mark values in bold if a significant difference between levels (95% CI)
tab_te3_hte8 <- sig_diff_func(tab_te3_hte8)

addtorow <- list(); addtorow$pos <- list(nrow(tab_te3_hte8)); addtorow$command <- c("")
ptable(tab_te3_hte8, file = "output/te3_hte8_cp_status.tex")

### 4.3.2 Table E4: Carbon tax status -------------------------------------
m_hte8a_eff <- amce(selected_profile ~ instrument + coverage + revenue_use + tax_status_factor + 
                      tax_status_factor:instrument + tax_status_factor:coverage + tax_status_factor:revenue_use, 
                    data = dat_conjoint %>% filter(outcome == "effective") %>%
                      mutate(tax_status_factor = factor(tax_status_factor,levels = c(0,1,2),labels = c("None","Under consideration","Adopted"))) %>% 
                      select(response_id, selected_profile, instrument, coverage, revenue_use, tax_status_factor) %>% na.omit(),
                    cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                    respondent.varying = "tax_status_factor")
m_hte8a_eff_sum <- summary(m_hte8a_eff)

m_hte8a_fea <- amce(selected_profile ~ instrument + coverage + revenue_use + tax_status_factor + 
                      tax_status_factor:instrument + tax_status_factor:coverage + tax_status_factor:revenue_use, 
                    data = dat_conjoint %>% filter(outcome == "feasible") %>%
                      mutate(tax_status_factor = factor(tax_status_factor,levels = c(0,1,2),labels = c("None","Under consideration","Adopted"))) %>% 
                      select(response_id, selected_profile, instrument, coverage, revenue_use, tax_status_factor) %>% na.omit(),
                    cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                    respondent.varying = "tax_status_factor")
m_hte8a_fea_sum <- summary(m_hte8a_fea)

m_hte8a_pro <- amce(selected_profile ~ instrument + coverage + revenue_use + tax_status_factor + 
                      tax_status_factor:instrument + tax_status_factor:coverage + tax_status_factor:revenue_use, 
                    data = dat_conjoint %>% filter(outcome == "propose") %>%
                      mutate(tax_status_factor = factor(tax_status_factor,levels = c(0,1,2),labels = c("None","Under consideration","Adopted"))) %>% 
                      select(response_id, selected_profile, instrument, coverage, revenue_use, tax_status_factor) %>% na.omit(),
                    cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                    respondent.varying = "tax_status_factor")
m_hte8a_pro_sum <- summary(m_hte8a_pro)

# Make table
mods_te4_hte8a <- list(m_hte8a_eff_sum,m_hte8a_fea_sum,m_hte8a_pro_sum)
mods_te4_hte8a_vecs <- tab_func(mods_te4_hte8a,amce_acie_resp = "resp")
mods_te4_hte8a_vecs$se_adj <- if_else(mods_te4_hte8a_vecs$se_adj=="(Inf)","\\multicolumn{1}{c}{\\textit{N.S.}}",mods_te4_hte8a_vecs$se_adj)
level_values_hte8a <- c("None","Under consideration","Adopted")

tab_te4_hte8a <- tab_hte_func(mods = mods_te4_hte8a, mods_vecs = mods_te4_hte8a_vecs, level_values = level_values_hte8a)

# Mark values in bold if a significant difference between levels (95% CI)
tab_te4_hte8a <- sig_diff_func(tab_te4_hte8a)

addtorow <- list(); addtorow$pos <- list(nrow(tab_te4_hte8a)); addtorow$command <- c("")
ptable(tab_te4_hte8a, file = "output/te4_hte8a_tax_status.tex")

### 4.3.3 Table E5: ETS status -------------------------------------
m_hte8b_eff <- amce(selected_profile ~ instrument + coverage + revenue_use + ets_status_factor + 
                      ets_status_factor:instrument + ets_status_factor:coverage + ets_status_factor:revenue_use, 
                    data = dat_conjoint %>% filter(outcome == "effective") %>%
                      mutate(ets_status_factor = factor(ets_status_factor,levels = c(0,1,2),labels = c("None","Under consideration","Adopted"))) %>% 
                      select(response_id, selected_profile, instrument, coverage, revenue_use, ets_status_factor) %>% na.omit(),
                    cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                    respondent.varying = "ets_status_factor")
m_hte8b_eff_sum <- summary(m_hte8b_eff)

m_hte8b_fea <- amce(selected_profile ~ instrument + coverage + revenue_use + ets_status_factor + 
                      ets_status_factor:instrument + ets_status_factor:coverage + ets_status_factor:revenue_use, 
                    data = dat_conjoint %>% filter(outcome == "feasible") %>%
                      mutate(ets_status_factor = factor(ets_status_factor,levels = c(0,1,2),labels = c("None","Under consideration","Adopted"))) %>% 
                      select(response_id, selected_profile, instrument, coverage, revenue_use, ets_status_factor) %>% na.omit(),
                    cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                    respondent.varying = "ets_status_factor")
m_hte8b_fea_sum <- summary(m_hte8b_fea)

m_hte8b_pro <- amce(selected_profile ~ instrument + coverage + revenue_use + ets_status_factor + 
                      ets_status_factor:instrument + ets_status_factor:coverage + ets_status_factor:revenue_use, 
                    data = dat_conjoint %>% filter(outcome == "propose") %>%
                      mutate(ets_status_factor = factor(ets_status_factor,levels = c(0,1,2),labels = c("None","Under consideration","Adopted"))) %>% 
                      select(response_id, selected_profile, instrument, coverage, revenue_use, ets_status_factor) %>% na.omit(),
                    cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                    respondent.varying = "ets_status_factor")
m_hte8b_pro_sum <- summary(m_hte8b_pro)

# Make table
mods_te5_hte8b <- list(m_hte8b_eff_sum,m_hte8b_fea_sum,m_hte8b_pro_sum)
mods_te5_hte8b_vecs <- tab_func(mods_te5_hte8b,amce_acie_resp = "resp")
mods_te5_hte8b_vecs$se_adj <- if_else(mods_te5_hte8b_vecs$se_adj=="(Inf)","\\multicolumn{1}{c}{\\textit{N.S.}}",mods_te5_hte8b_vecs$se_adj)
level_values_hte8b <- c("None","Under consideration","Adopted")

tab_te5_hte8b <- tab_hte_func(mods = mods_te5_hte8b, mods_vecs = mods_te5_hte8b_vecs, level_values = level_values_hte8b)

# Mark values in bold if a significant difference between levels (95% CI)
tab_te5_hte8b <- sig_diff_func(tab_te5_hte8b)

addtorow <- list(); addtorow$pos <- list(nrow(tab_te5_hte8b)); addtorow$command <- c("")
ptable(tab_te5_hte8b, file = "output/te5_hte8b_ets_status.tex")

## 4.4 Demand for carbon pricing -------------------------------------

### 4.4.1 Table E6: Pressure from domestic vs. international sources -------------------------------------
m_hte2_eff <- amce(selected_profile ~ instrument + coverage + revenue_use + driver_domestic_international + 
                     driver_domestic_international:instrument + driver_domestic_international:coverage + driver_domestic_international:revenue_use, 
                   data = dat_conjoint %>% filter(outcome == "effective") %>% 
                     select(response_id, selected_profile, instrument, coverage, revenue_use, driver_domestic_international) %>% na.omit() %>%
                     mutate(driver_domestic_international = factor(driver_domestic_international,levels = c(-1,0,1),
                                                                   labels = c("Mostly domestic",
                                                                              "Both international and domestic equally",
                                                                              "Mostly international"))),
                   cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                   respondent.varying = "driver_domestic_international")
m_hte2_eff_sum <- summary(m_hte2_eff)

m_hte2_fea <- amce(selected_profile ~ instrument + coverage + revenue_use + driver_domestic_international + 
                     driver_domestic_international:instrument + driver_domestic_international:coverage + driver_domestic_international:revenue_use, 
                   data = dat_conjoint %>% filter(outcome == "feasible") %>% 
                     select(response_id, selected_profile, instrument, coverage, revenue_use, driver_domestic_international) %>% na.omit() %>%
                     mutate(driver_domestic_international = factor(driver_domestic_international,levels = c(-1,0,1),
                                                                   labels = c("Mostly domestic",
                                                                              "Both international and domestic equally",
                                                                              "Mostly international"))),
                   cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                   respondent.varying = "driver_domestic_international")
m_hte2_fea_sum <- summary(m_hte2_fea)


m_hte2_pro <- amce(selected_profile ~ instrument + coverage + revenue_use + driver_domestic_international + 
                     driver_domestic_international:instrument + driver_domestic_international:coverage + driver_domestic_international:revenue_use, 
                   data = dat_conjoint %>% filter(outcome == "propose") %>% 
                     select(response_id, selected_profile, instrument, coverage, revenue_use, driver_domestic_international) %>% na.omit() %>%
                     mutate(driver_domestic_international = factor(driver_domestic_international,levels = c(-1,0,1),
                                                                   labels = c("Mostly domestic",
                                                                              "Both international and domestic equally",
                                                                              "Mostly international"))),
                   cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                   respondent.varying = "driver_domestic_international")
m_hte2_pro_sum <- summary(m_hte2_pro)

mods_te6_hte2 <- list(m_hte2_eff_sum,m_hte2_fea_sum,m_hte2_pro_sum)
mods_te6_hte2_vecs <- tab_func(mods_te6_hte2,amce_acie_resp = "resp")
mods_te6_hte2_vecs$se_adj <- if_else(mods_te6_hte2_vecs$se_adj=="(Inf)","\\multicolumn{1}{c}{\\textit{N.S.}}",mods_te6_hte2_vecs$se_adj)
level_values_hte2 <- c("Domestic","Both equally","International")

tab_te6_hte2 <- tab_hte_func(mods = mods_te6_hte2, mods_vecs = mods_te6_hte2_vecs, level_values = level_values_hte2)

# Mark values in bold if a significant difference between levels (95% CI)
tab_te6_hte2 <- sig_diff_func(tab_te6_hte2)

addtorow <- list(); addtorow$pos <- list(nrow(tab_te6_hte2)); addtorow$command <- c("")
ptable(tab_te6_hte2, file = "output/te6_hte2_domestic_international.tex")

### 4.4.2 Table E7: Pressure from public vs. private sector -------------------------------------
m_hte3_eff <- amce(selected_profile ~ instrument + coverage + revenue_use + driver_private_public + 
                     driver_private_public:instrument + driver_private_public:coverage + driver_private_public:revenue_use, 
                   data = dat_conjoint %>% filter(outcome == "effective") %>% 
                     select(response_id, selected_profile, instrument, coverage, revenue_use, driver_private_public) %>% na.omit() %>%
                     mutate(driver_private_public = factor(driver_private_public,levels = c(-1,0,1),
                                                           labels = c("Mostly the public",
                                                                      "Both the private sector and the public equally",
                                                                      "Mostly the private sector"))),
                   cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                   respondent.varying = "driver_private_public")
m_hte3_eff_sum <- summary(m_hte3_eff)

m_hte3_fea <- amce(selected_profile ~ instrument + coverage + revenue_use + driver_private_public + 
                     driver_private_public:instrument + driver_private_public:coverage + driver_private_public:revenue_use, 
                   data = dat_conjoint %>% filter(outcome == "feasible") %>% 
                     select(response_id, selected_profile, instrument, coverage, revenue_use, driver_private_public) %>% na.omit() %>%
                     mutate(driver_private_public = factor(driver_private_public,levels = c(-1,0,1),
                                                           labels = c("Mostly the public",
                                                                      "Both the private sector and the public equally",
                                                                      "Mostly the private sector"))),
                   cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                   respondent.varying = "driver_private_public")
m_hte3_fea_sum <- summary(m_hte3_fea)


m_hte3_pro <- amce(selected_profile ~ instrument + coverage + revenue_use + driver_private_public + 
                     driver_private_public:instrument + driver_private_public:coverage + driver_private_public:revenue_use, 
                   data = dat_conjoint %>% filter(outcome == "propose") %>% 
                     select(response_id, selected_profile, instrument, coverage, revenue_use, driver_private_public) %>% na.omit() %>%
                     mutate(driver_private_public = factor(driver_private_public,levels = c(-1,0,1),
                                                           labels = c("Mostly the public",
                                                                      "Both the private sector and the public equally",
                                                                      "Mostly the private sector"))),
                   cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                   respondent.varying = "driver_private_public")
m_hte3_pro_sum <- summary(m_hte3_pro)

mods_te7_hte3 <- list(m_hte3_eff_sum,m_hte3_fea_sum,m_hte3_pro_sum)
mods_te7_hte3_vecs <- tab_func(mods_te7_hte3,amce_acie_resp = "resp")
mods_te7_hte3_vecs$se_adj <- if_else(mods_te7_hte3_vecs$se_adj=="(Inf)","\\multicolumn{1}{c}{\\textit{N.S.}}",mods_te7_hte3_vecs$se_adj)
level_values_hte3 <- c("Public","Both equally","Private")

tab_te7_hte3 <- tab_hte_func(mods = mods_te7_hte3, mods_vecs = mods_te7_hte3_vecs, level_values = level_values_hte3)

# Mark values in bold if a significant difference between levels (95% CI)
tab_te7_hte3 <- sig_diff_func(tab_te7_hte3)

addtorow <- list(); addtorow$pos <- list(nrow(tab_te7_hte3)); addtorow$command <- c("")
ptable(tab_te7_hte3, file = "output/te7_hte3_private_public.tex")

### 4.4.3 Table E8: Political priority on adaptation vs. mitigation -------------------------------------
m_hte4_eff <- amce(selected_profile ~ instrument + coverage + revenue_use + adaptation_mitigation + 
                     adaptation_mitigation:instrument + adaptation_mitigation:coverage + adaptation_mitigation:revenue_use, 
                   data = dat_conjoint %>% filter(outcome == "effective") %>% 
                     select(response_id, selected_profile, instrument, coverage, revenue_use, adaptation_mitigation) %>% na.omit() %>%
                     mutate(adaptation_mitigation = factor(adaptation_mitigation,levels = c(-1,0,1),
                                                           labels = c("Mostly adaptation",
                                                                      "Both mitigation and adaptation equally",
                                                                      "Mostly mitigation"))),
                   cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                   respondent.varying = "adaptation_mitigation")
m_hte4_eff_sum <- summary(m_hte4_eff)

m_hte4_fea <- amce(selected_profile ~ instrument + coverage + revenue_use + adaptation_mitigation + 
                     adaptation_mitigation:instrument + adaptation_mitigation:coverage + adaptation_mitigation:revenue_use, 
                   data = dat_conjoint %>% filter(outcome == "feasible") %>% 
                     select(response_id, selected_profile, instrument, coverage, revenue_use, adaptation_mitigation) %>% na.omit() %>%
                     mutate(adaptation_mitigation = factor(adaptation_mitigation,levels = c(-1,0,1),
                                                           labels = c("Mostly adaptation",
                                                                      "Both mitigation and adaptation equally",
                                                                      "Mostly mitigation"))),
                   cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                   respondent.varying = "adaptation_mitigation")
m_hte4_fea_sum <- summary(m_hte4_fea)


m_hte4_pro <- amce(selected_profile ~ instrument + coverage + revenue_use + adaptation_mitigation + 
                     adaptation_mitigation:instrument + adaptation_mitigation:coverage + adaptation_mitigation:revenue_use, 
                   data = dat_conjoint %>% filter(outcome == "propose") %>% 
                     select(response_id, selected_profile, instrument, coverage, revenue_use, adaptation_mitigation) %>% na.omit() %>%
                     mutate(adaptation_mitigation = factor(adaptation_mitigation,levels = c(-1,0,1),
                                                           labels = c("Mostly adaptation",
                                                                      "Both mitigation and adaptation equally",
                                                                      "Mostly mitigation"))),
                   cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                   respondent.varying = "adaptation_mitigation")
m_hte4_pro_sum <- summary(m_hte4_pro)

mods_te8_hte4 <- list(m_hte4_eff_sum,m_hte4_fea_sum,m_hte4_pro_sum)
mods_te8_hte4_vecs <- tab_func(mods_te8_hte4,amce_acie_resp = "resp")
mods_te8_hte4_vecs$se_adj <- if_else(mods_te8_hte4_vecs$se_adj=="(Inf)","\\multicolumn{1}{c}{\\textit{N.S.}}",mods_te8_hte4_vecs$se_adj)
level_values_hte4 <- c("Adapatation","Both equally","Mitigation")

tab_te8_hte4 <- tab_hte_func(mods = mods_te8_hte4, mods_vecs = mods_te8_hte4_vecs, level_values = level_values_hte4)

# Mark values in bold if a significant difference between levels (95% CI)
tab_te8_hte4 <- sig_diff_func(tab_te8_hte4)

addtorow <- list(); addtorow$pos <- list(nrow(tab_te8_hte4)); addtorow$command <- c("")
ptable(tab_te8_hte4, file = "output/te8_hte4_adaptation_mitigation.tex")

## 4.5 Country context -------------------------------------

### 4.5.1 Table E9: Electoral democracy -------------------------------------
m_hte5_eff <- amce(selected_profile ~ instrument + coverage + revenue_use + v2x_polyarchy + 
                     v2x_polyarchy:instrument + v2x_polyarchy:coverage + v2x_polyarchy:revenue_use, 
                   data = dat_conjoint %>% filter(outcome == "effective") %>% 
                     select(response_id, selected_profile, instrument, coverage, revenue_use, v2x_polyarchy) %>% na.omit(),
                   cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                   respondent.varying = "v2x_polyarchy")
m_hte5_eff_sum <- summary(m_hte5_eff,covariate.values = list("v2x_polyarchy" = quantile(m_hte5_eff$data$v2xpolyarchy,probs = c(0.2,0.8))))

m_hte5_fea <- amce(selected_profile ~ instrument + coverage + revenue_use + v2x_polyarchy + 
                     v2x_polyarchy:instrument + v2x_polyarchy:coverage + v2x_polyarchy:revenue_use, 
                   data = dat_conjoint %>% filter(outcome == "feasible") %>% 
                     select(response_id, selected_profile, instrument, coverage, revenue_use, v2x_polyarchy) %>% na.omit(),
                   cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                   respondent.varying = "v2x_polyarchy")
m_hte5_fea_sum <- summary(m_hte5_fea,covariate.values = list("v2x_polyarchy" = quantile(m_hte5_fea$data$v2xpolyarchy,probs = c(0.2,0.8))))


m_hte5_pro <- amce(selected_profile ~ instrument + coverage + revenue_use + v2x_polyarchy + 
                     v2x_polyarchy:instrument + v2x_polyarchy:coverage + v2x_polyarchy:revenue_use, 
                   data = dat_conjoint %>% filter(outcome == "propose") %>% 
                     select(response_id, selected_profile, instrument, coverage, revenue_use, v2x_polyarchy) %>% na.omit(),
                   cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                   respondent.varying = "v2x_polyarchy")
m_hte5_pro_sum <- summary(m_hte5_pro,covariate.values = list("v2x_polyarchy" = quantile(m_hte5_pro$data$v2xpolyarchy,probs = c(0.2,0.8))))

mods_te9_hte5 <- list(m_hte5_eff_sum,m_hte5_fea_sum,m_hte5_pro_sum)
mods_te9_hte5_vecs <- tab_func(mods_te9_hte5,amce_acie_resp = "resp")
mods_te9_hte5_vecs$se_adj <- if_else(mods_te9_hte5_vecs$se_adj=="(Inf)","\\multicolumn{1}{c}{\\textit{N.S.}}",mods_te9_hte5_vecs$se_adj)
level_values_hte5 <- paste0("P\\_{",str_remove(mods_te9_hte5[[1]]$table_values_amce[,"Level.Value"],pattern = "%"),"\\%}")

tab_te9_hte5 <- tab_hte_func(mods = mods_te9_hte5, mods_vecs = mods_te9_hte5_vecs, level_values = level_values_hte5)

# Mark values in bold if a significant difference between levels (95% CI)
tab_te9_hte5 <- sig_diff_func(tab_te9_hte5)

addtorow <- list(); addtorow$pos <- list(nrow(tab_te9_hte5)); addtorow$command <- c("")
ptable(tab_te9_hte5, file = "output/te9_hte5_democracy.tex")

### 4.5.2 Table E10: GDP per capita -------------------------------------
m_hte6_eff <- amce(selected_profile ~ instrument + coverage + revenue_use + gdp_pc + 
                     gdp_pc:instrument + gdp_pc:coverage + gdp_pc:revenue_use, 
                   data = dat_conjoint %>% filter(outcome == "effective") %>% 
                     select(response_id, selected_profile, instrument, coverage, revenue_use, gdp_pc) %>% na.omit(),
                   cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                   respondent.varying = "gdp_pc")
m_hte6_eff_sum <- summary(m_hte6_eff,covariate.values = list("gdp_pc" = quantile(m_hte6_eff$data$gdppc,probs = c(0.2,0.8))))

m_hte6_fea <- amce(selected_profile ~ instrument + coverage + revenue_use + gdp_pc + 
                     gdp_pc:instrument + gdp_pc:coverage + gdp_pc:revenue_use, 
                   data = dat_conjoint %>% filter(outcome == "feasible") %>% 
                     select(response_id, selected_profile, instrument, coverage, revenue_use, gdp_pc) %>% na.omit(),
                   cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                   respondent.varying = "gdp_pc")
m_hte6_fea_sum <- summary(m_hte6_fea,covariate.values = list("gdp_pc" = quantile(m_hte6_fea$data$gdppc,probs = c(0.2,0.8))))


m_hte6_pro <- amce(selected_profile ~ instrument + coverage + revenue_use + gdp_pc + 
                     gdp_pc:instrument + gdp_pc:coverage + gdp_pc:revenue_use, 
                   data = dat_conjoint %>% filter(outcome == "propose") %>% 
                     select(response_id, selected_profile, instrument, coverage, revenue_use, gdp_pc) %>% na.omit(),
                   cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                   respondent.varying = "gdp_pc")
m_hte6_pro_sum <- summary(m_hte6_pro,covariate.values = list("gdp_pc" = quantile(m_hte6_pro$data$gdppc,probs = c(0.2,0.8))))

mods_te10_hte6 <- list(m_hte6_eff_sum,m_hte6_fea_sum,m_hte6_pro_sum)
mods_te10_hte6_vecs <- tab_func(mods_te10_hte6,amce_acie_resp = "resp")
mods_te10_hte6_vecs$se_adj <- if_else(mods_te10_hte6_vecs$se_adj=="(Inf)","\\multicolumn{1}{c}{\\textit{N.S.}}",mods_te10_hte6_vecs$se_adj)
level_values_hte6 <- paste0("P\\_{",str_remove(mods_te10_hte6[[1]]$table_values_amce[,"Level.Value"],pattern = "%"),"\\%}")

tab_te10_hte6 <- tab_hte_func(mods = mods_te10_hte6, mods_vecs = mods_te10_hte6_vecs, level_values = level_values_hte6)

# Mark values in bold if a significant difference between levels (95% CI)
tab_te10_hte6 <- sig_diff_func(tab_te10_hte6)

addtorow <- list(); addtorow$pos <- list(nrow(tab_te10_hte6)); addtorow$command <- c("")
ptable(tab_te10_hte6, file = "output/te10_hte6_gdp_pc.tex")

### 4.5.3 Table E11: ODA per capita -------------------------------------
m_hte7_eff <- amce(selected_profile ~ instrument + coverage + revenue_use + oda_pc + 
                     oda_pc:instrument + oda_pc:coverage + oda_pc:revenue_use, 
                   data = dat_conjoint %>% filter(outcome == "effective") %>% 
                     select(response_id, selected_profile, instrument, coverage, revenue_use, oda_pc) %>% na.omit(),
                   cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                   respondent.varying = "oda_pc")
m_hte7_eff_sum <- summary(m_hte7_eff,covariate.values = list("oda_pc" = quantile(m_hte7_eff$data$odapc,probs = c(0.2,0.8))))

m_hte7_fea <- amce(selected_profile ~ instrument + coverage + revenue_use + oda_pc + 
                     oda_pc:instrument + oda_pc:coverage + oda_pc:revenue_use, 
                   data = dat_conjoint %>% filter(outcome == "feasible") %>% 
                     select(response_id, selected_profile, instrument, coverage, revenue_use, oda_pc) %>% na.omit(),
                   cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                   respondent.varying = "oda_pc")
m_hte7_fea_sum <- summary(m_hte7_fea,covariate.values = list("oda_pc" = quantile(m_hte7_fea$data$odapc,probs = c(0.2,0.8))))

m_hte7_pro <- amce(selected_profile ~ instrument + coverage + revenue_use + oda_pc + 
                     oda_pc:instrument + oda_pc:coverage + oda_pc:revenue_use, 
                   data = dat_conjoint %>% filter(outcome == "propose") %>% 
                     select(response_id, selected_profile, instrument, coverage, revenue_use, oda_pc) %>% na.omit(),
                   cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                   respondent.varying = "oda_pc")
m_hte7_pro_sum <- summary(m_hte7_pro,covariate.values = list("oda_pc" = quantile(m_hte7_pro$data$odapc,probs = c(0.2,0.8))))

mods_te11_hte7 <- list(m_hte7_eff_sum,m_hte7_fea_sum,m_hte7_pro_sum)
mods_te11_hte7_vecs <- tab_func(mods_te11_hte7,amce_acie_resp = "resp")
mods_te11_hte7_vecs$se_adj <- if_else(mods_te11_hte7_vecs$se_adj=="(Inf)","\\multicolumn{1}{c}{\\textit{N.S.}}",mods_te11_hte7_vecs$se_adj)
level_values_hte7 <- paste0("P\\_{",str_remove(mods_te11_hte7[[1]]$table_values_amce[,"Level.Value"],pattern = "%"),"\\%}")

tab_te11_hte7 <- tab_hte_func(mods = mods_te11_hte7, mods_vecs = mods_te11_hte7_vecs, level_values = level_values_hte7)

# Mark values in bold if a significant difference between levels (95% CI)
tab_te11_hte7 <- sig_diff_func(tab_te11_hte7)

addtorow <- list(); addtorow$pos <- list(nrow(tab_te11_hte7)); addtorow$command <- c("")
ptable(tab_te11_hte7, file = "output/te11_hte7_oda_pc.tex")

## 4.6 Respondent characteristics -------------------------------------

### 4.6.1 Table E12: Respondent carbon pricing support ------------------------------------
m_hte1_eff <- amce(selected_profile ~ instrument + coverage + revenue_use + cp_support + 
                     cp_support:instrument + cp_support:coverage + cp_support:revenue_use, 
                   data = dat_conjoint %>% filter(outcome == "effective") %>% 
                     select(response_id, selected_profile, instrument, coverage, revenue_use, cp_support) %>% na.omit(),
                   cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                   respondent.varying = "cp_support")
m_hte1_eff_sum <- summary(m_hte1_eff,covariate.values = list("cp_support" = quantile(m_hte1_eff$data$cpsupport,probs = c(0.2,0.8))))

m_hte1_fea <- amce(selected_profile ~ instrument + coverage + revenue_use + cp_support + 
                     cp_support:instrument + cp_support:coverage + cp_support:revenue_use, 
                   data = dat_conjoint %>% filter(outcome == "feasible") %>% 
                     select(response_id, selected_profile, instrument, coverage, revenue_use, cp_support) %>% na.omit(),
                   cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                   respondent.varying = "cp_support")
m_hte1_fea_sum <- summary(m_hte1_fea,covariate.values = list("cp_support" = quantile(m_hte1_fea$data$cpsupport,probs = c(0.2,0.8))))


m_hte1_pro <- amce(selected_profile ~ instrument + coverage + revenue_use + cp_support + 
                     cp_support:instrument + cp_support:coverage + cp_support:revenue_use, 
                   data = dat_conjoint %>% filter(outcome == "propose") %>% 
                     select(response_id, selected_profile, instrument, coverage, revenue_use, cp_support) %>% na.omit(),
                   cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                   respondent.varying = "cp_support")
m_hte1_pro_sum <- summary(m_hte1_pro,covariate.values = list("cp_support" = quantile(m_hte1_pro$data$cpsupport,probs = c(0.2,0.8))))

mods_te12_hte1 <- list(m_hte1_eff_sum,m_hte1_fea_sum,m_hte1_pro_sum)
mods_te12_hte1_vecs <- tab_func(mods_te12_hte1,amce_acie_resp = "resp")
mods_te12_hte1_vecs$se_adj <- if_else(mods_te12_hte1_vecs$se_adj=="(Inf)","\\multicolumn{1}{c}{\\textit{N.S.}}",mods_te12_hte1_vecs$se_adj)
level_values_hte1 <- paste0("P\\_{",str_remove(mods_te12_hte1[[1]]$table_values_amce[,"Level.Value"],pattern = "%"),"\\%}")

tab_te12_hte1 <- tab_hte_func(mods = mods_te12_hte1, mods_vecs = mods_te12_hte1_vecs, level_values = level_values_hte1)

# Mark values in bold if a significant difference between levels (95% CI)
tab_te12_hte1 <- sig_diff_func(tab_te12_hte1)

addtorow <- list(); addtorow$pos <- list(nrow(tab_te12_hte1)); addtorow$command <- c("")
ptable(tab_te12_hte1, file = "output/te12_hte1_cp_support.tex")

### 4.6.2 Table E13: International vs. domestic respondent ------------------------------------
m_hte9_eff <- amce(selected_profile ~ instrument + coverage + revenue_use + location + 
                     location:instrument + location:coverage + location:revenue_use, 
                   data = dat_conjoint %>% filter(outcome == "effective") %>%
                     select(response_id, selected_profile, instrument, coverage, revenue_use, location) %>% na.omit() %>%
                     mutate(location = factor(location,levels = c("Outside country","Inside country"))),
                   cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                   respondent.varying = "location")
m_hte9_eff_sum <- summary(m_hte9_eff)

m_hte9_fea <- amce(selected_profile ~ instrument + coverage + revenue_use + location + 
                     location:instrument + location:coverage + location:revenue_use, 
                   data = dat_conjoint %>% filter(outcome == "feasible") %>%
                     select(response_id, selected_profile, instrument, coverage, revenue_use, location) %>% na.omit() %>%
                     mutate(location = factor(location,levels = c("Outside country","Inside country"))),
                   cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                   respondent.varying = "location")
m_hte9_fea_sum <- summary(m_hte9_fea)


m_hte9_pro <- amce(selected_profile ~ instrument + coverage + revenue_use + location + 
                     location:instrument + location:coverage + location:revenue_use, 
                   data = dat_conjoint %>% filter(outcome == "propose") %>%
                     select(response_id, selected_profile, instrument, coverage, revenue_use, location) %>% na.omit() %>%
                     mutate(location = factor(location,levels = c("Outside country","Inside country"))),
                   cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                   respondent.varying = "location")
m_hte9_pro_sum <- summary(m_hte9_pro)

mods_te13_hte9 <- list(m_hte9_eff_sum,m_hte9_fea_sum,m_hte9_pro_sum)
mods_te13_hte9_vecs <- tab_func(mods_te13_hte9,amce_acie_resp = "resp")
mods_te13_hte9_vecs$se_adj <- if_else(mods_te13_hte9_vecs$se_adj=="(Inf)","\\multicolumn{1}{c}{\\textit{N.S.}}",mods_te13_hte9_vecs$se_adj)
level_values_hte9 <- c("Outside country","Inside country")

tab_te13_hte9 <- tab_hte_func(mods = mods_te13_hte9, mods_vecs = mods_te13_hte9_vecs, level_values = level_values_hte9)

# Mark values in bold if a significant difference between levels (95% CI)
tab_te13_hte9 <- sig_diff_func(tab_te13_hte9)

addtorow <- list(); addtorow$pos <- list(nrow(tab_te13_hte9)); addtorow$command <- c("")
ptable(tab_te13_hte9, file = "output/te13_hte9_location.tex")

### 4.6.3 Table E14: Government vs. non-government position ------------------------------------
m_rq6_eff <- amce(selected_profile ~ instrument + coverage + revenue_use + in_government + 
                    in_government:instrument + in_government:coverage + in_government:revenue_use, 
                  data = dat_conjoint %>% filter(outcome == "effective") %>% mutate(in_government = factor(in_government)), cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                  respondent.varying = "in_government")
m_rq6_eff_sum <- summary(m_rq6_eff)

m_rq6_fea <- amce(selected_profile ~ instrument + coverage + revenue_use + in_government + 
                    in_government:instrument + in_government:coverage + in_government:revenue_use, 
                  data = dat_conjoint %>% filter(outcome == "feasible") %>% mutate(in_government = factor(in_government)), cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                  respondent.varying = "in_government")
m_rq6_fea_sum <- summary(m_rq6_fea)

m_rq6_pro <- amce(selected_profile ~ instrument + coverage + revenue_use + in_government + 
                    in_government:instrument + in_government:coverage + in_government:revenue_use, 
                  data = dat_conjoint %>% filter(outcome == "propose") %>% mutate(in_government = factor(in_government)), cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                  respondent.varying = "in_government")
m_rq6_pro_sum <- summary(m_rq6_pro)

mods_te14_rq6 <- list(m_rq6_eff_sum,m_rq6_fea_sum,m_rq6_pro_sum)
mods_te14_rq6_vecs <- tab_func(mods_te14_rq6,amce_acie_resp = "resp")
mods_te14_rq6_vecs$se_adj <- if_else(mods_te14_rq6_vecs$se_adj=="(Inf)","\\multicolumn{1}{c}{\\textit{N.S.}}",mods_te14_rq6_vecs$se_adj)
level_values_rq6 <- c("Non-government","Government")

tab_te14_rq6 <- tab_hte_func(mods = mods_te14_rq6, mods_vecs = mods_te14_rq6_vecs, level_values = level_values_rq6)

# Mark values in bold if a significant difference between levels (95% CI)
tab_te14_rq6 <- sig_diff_func(tab_te14_rq6)

addtorow <- list(); addtorow$pos <- list(nrow(tab_te14_rq6)); addtorow$command <- c("")
ptable(tab_te14_rq6, file = "output/te14_rq6_position.tex")

### 4.6.4 Table E15: Career stage ------------------------------------
m_hte11_eff <- amce(selected_profile ~ instrument + coverage + revenue_use + years_experience + 
                      years_experience:instrument + years_experience:coverage + years_experience:revenue_use, 
                    data = dat_conjoint %>% filter(outcome == "effective") %>%
                      mutate(years_experience = factor(years_experience)) %>%
                      select(response_id, selected_profile, instrument, coverage, revenue_use, years_experience) %>% na.omit(),
                    cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                    respondent.varying = "years_experience")
m_hte11_eff_sum <- summary(m_hte11_eff)

m_hte11_fea <- amce(selected_profile ~ instrument + coverage + revenue_use + years_experience + 
                      years_experience:instrument + years_experience:coverage + years_experience:revenue_use, 
                    data = dat_conjoint %>% filter(outcome == "feasible") %>%
                      mutate(years_experience = factor(years_experience)) %>%
                      select(response_id, selected_profile, instrument, coverage, revenue_use, years_experience) %>% na.omit(),
                    cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                    respondent.varying = "years_experience")
m_hte11_fea_sum <- summary(m_hte11_fea)


m_hte11_pro <- amce(selected_profile ~ instrument + coverage + revenue_use + years_experience + 
                      years_experience:instrument + years_experience:coverage + years_experience:revenue_use, 
                    data = dat_conjoint %>% filter(outcome == "propose") %>%
                      mutate(years_experience = factor(years_experience)) %>%
                      select(response_id, selected_profile, instrument, coverage, revenue_use, years_experience) %>% na.omit(),
                    cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                    respondent.varying = "yearsexperience")
m_hte11_pro_sum <- summary(m_hte11_pro)

mods_te15_hte11 <- list(m_hte11_eff_sum,m_hte11_fea_sum,m_hte11_pro_sum)
mods_te15_hte11_vecs <- tab_func(mods_te15_hte11,amce_acie_resp = "resp")
mods_te15_hte11_vecs$se_adj <- if_else(mods_te15_hte11_vecs$se_adj=="(Inf)","\\multicolumn{1}{c}{\\textit{N.S.}}",mods_te15_hte11_vecs$se_adj)
level_values_hte11 <- c("1--10 years","11--20 years","21--30 years","More than 30 years")

tab_te15_hte11 <- tab_hte_func(mods = mods_te15_hte11, mods_vecs = mods_te15_hte11_vecs, level_values = level_values_hte11)

# Mark values in bold if a significant difference between levels (95% CI)
tab_te15_hte11 <- sig_diff_func(tab_te15_hte11)

addtorow <- list(); addtorow$pos <- list(nrow(tab_te15_hte11)); addtorow$command <- c("")
ptable(tab_te15_hte11, file = "output/te15_hte11_career_stage.tex")

### 4.6.4 Table E16: Gender ------------------------------------
m_hte12_eff <- amce(selected_profile ~ instrument + coverage + revenue_use + gender + 
                      gender:instrument + gender:coverage + gender:revenue_use, 
                    data = dat_conjoint %>% filter(outcome == "effective") %>%
                      mutate(gender = factor(gender, levels = c("Male","Female"))) %>%
                      select(response_id, selected_profile, instrument, coverage, revenue_use, gender) %>% na.omit(),
                    cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                    respondent.varying = "gender")
m_hte12_eff_sum <- summary(m_hte12_eff)

m_hte12_fea <- amce(selected_profile ~ instrument + coverage + revenue_use + gender + 
                      gender:instrument + gender:coverage + gender:revenue_use, 
                    data = dat_conjoint %>% filter(outcome == "feasible") %>%
                      mutate(gender = factor(gender, levels = c("Male","Female"))) %>%
                      select(response_id, selected_profile, instrument, coverage, revenue_use, gender) %>% na.omit(),
                    cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                    respondent.varying = "gender")
m_hte12_fea_sum <- summary(m_hte12_fea)


m_hte12_pro <- amce(selected_profile ~ instrument + coverage + revenue_use + gender + 
                      gender:instrument + gender:coverage + gender:revenue_use, 
                    data = dat_conjoint %>% filter(outcome == "propose") %>%
                      mutate(gender = factor(gender, levels = c("Male","Female"))) %>%
                      select(response_id, selected_profile, instrument, coverage, revenue_use, gender) %>% na.omit(),
                    cluster = TRUE, respondent.id = "response_id",design = conjointdesign, 
                    respondent.varying = "gender")
m_hte12_pro_sum <- summary(m_hte12_pro)

mods_te16_hte12 <- list(m_hte12_eff_sum,m_hte12_fea_sum,m_hte12_pro_sum)
mods_te16_hte12_vecs <- tab_func(mods_te16_hte12,amce_acie_resp = "resp")
mods_te16_hte12_vecs$se_adj <- if_else(mods_te16_hte12_vecs$se_adj=="(Inf)","\\multicolumn{1}{c}{\\textit{N.S.}}",mods_te16_hte12_vecs$se_adj)
level_values_hte12 <- c("Male","Female")

tab_te16_hte12 <- tab_hte_func(mods = mods_te16_hte12, mods_vecs = mods_te16_hte12_vecs, level_values = level_values_hte12)

# Mark values in bold if a significant difference between levels (95% CI)
tab_te16_hte12 <- sig_diff_func(tab_te16_hte12)

addtorow <- list(); addtorow$pos <- list(nrow(tab_te16_hte12)); addtorow$command <- c("")
ptable(tab_te16_hte12, file = "output/te16_hte12_gender.tex")

### 4.6.5 Fossil fuel sector respondent -----------------------------
# Pre-registered but not run due to low number of fossil fuel sector respondents